Quantum Mechanical Breakdown of Perfect Homogeneity in Reheating After Inflation 
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In the context of quantum fields in time dependent classical backgrounds, we notice that the 
number of created particles with a given momentum largely deviates about its mean value. Guided 
with this observation we use a complete orthonormal family of localized wave packets to calculate 
the deviations in the number and energy densities of particles produced in a volume of a given size 
during reheating. It turns out that at the end of reheating there exists (in general tiny) spatial 
variations in these densities on Hubble length scales over which local interactions are incapable of 
restoring homogeneity. This signals the destruction of perfect homogeneity attained after inflation 
due to the quantum nature of particle production process in reheating. 



I. INTRODUCTION 

In scalar field driven inflationary models, the universe 
expands (almost) exponentially due to the scalar poten- 
tial energy density acting as an effective cosmological 
constant. As a result, a small pre-inflationary causal 
patch enlarges to encompass the whole observed universe. 
With this notion, inflation solves many puzzles of the 
standard cosmological model like homogeneity, isotropy 
and monopolc problems. 

As inflation ends, one finds an extremely smooth, flat 
universe which is practically at zero temperature and 
filled only with a coherently oscillating inflaton field of 
decreasing amplitude about the minimum of its poten- 
tial. Almost all matter populating the universe in the 
subsequent radiation era should have been produced by 
the decay of the inflaton field. These decay products later 
become thermalized by collisions and further decays and 
this marks the beginning of the usual hot stage. 

To have a complete cosmic history of the universe, it is 
crucial to understand reheating after inflation in detail. 
The elementary theory of reheating (for a review see e.g. 
[1]) is based on perturbation theory which has obvious 
limitations. More recently, the importance of the para- 
metric resonance effects on the decay of the inflaton field 
is recognized in [2-10] (for a recent review see, e.g. [11]). 
Especially it is shown in [8] that in many models reheat- 
ing actually starts by a decay due to broad parametric 
resonance, called preheating. The particle creation ef- 
fects in the broad parametric resonance regime is stud- 
ied in [10], which shows that the whole process becomes 
stochastic due to the expansion of the universe but one 
still gets exponentially growing occupation numbers. 

Consider the universe just after the inflation at the 
beginning of the reheating stage, which has already ac- 
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quired a huge size. Viewing this moment as the starting 
point of the subsequent evolution, causality requires that 
the particle creation and thermalization processes occur 
independently in each succeeding Hubble volume, which 
(ignoring the slight difference between the horizon and 
the Hubble distances) form causally disjoint regions dur- 
ing reheating. Naively thinking, this does not imply a 
breakdown in homogeneity since in each region we have 
identical copies of the same oscillating scalar field. How- 
ever, since the decay process is quantum mechanical in 
nature, the number of created particles in each volume 
is expected to fluctuate about a mean value which may 
induce (in general small) density inhomogeneities. The 
aim of this work is to calculate these variations by analyz- 
ing the decay process of the inflaton field using localized 
wave packets. 

The plan of the paper is as follows. In the follow- 
ing section, we review the well known properties of the 
quantum mechanical harmonic oscillator with a time de- 
pendent frequency. We point out some salient features of 
this system which will be important in the following sec- 
tions. In section III, we study the general framework of 
particle production by time dependent external classical 
fields using a complete orthonormal family of localized 
wave packets and determine the deviations in the num- 
ber and the energy densities of particles produced in a 
given volume. In section IV , we apply these results to a 
generic inflationary scenario and estimate the amount of 
inhomogeneities. We conclude with a brief review of our 
findings in V. 



II. HARMONIC OSCILLATOR WITH TIME 
DEPENDENT FREQUENCY 

In this section we review the basic properties of the 
quantum mechanical harmonic oscillator with time de- 
pendent frequency lo. As it is well known, this sys- 
tem mimics the essential features of the particle creation 
process during reheating. For a unit mass particle, the 
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Hamiltonian of the oscillator is given by 
1 



(1) 



where q and p denote the position and the momentum 
operators obeying 

[q,p]=i- ( 2 ) 
In the Hciscnberg picture these operators obey 

q = p, p = -uq, (3) 

where dot denotes time derivative. To solve the system 
one may introduce time dependent ladder operators a 
and a* via 



q = 



ij LO(lt 



2lj 



at 



i J udt 



2lj 



(4) 



In taking time derivative of q, one wishes to treat a, at 
and uj as if they were time independent quantities (like 
in the usual harmonic oscillator) so that 



P : 



-if udt 



i J udt 



'2uj V2uj 
This can be achieved by imposing 



(5) 



u 2i f' udt f 

a = — e J a\ 
2u> 



at 



2uj 



-2ifu,dt a (6) 



The nice feature of defining the ladder operators in this 
fashion is that the canonical commutation relation (2) 
gives 



[a, a f ] = 1, 
and the Hamiltonian becomes 



H 



t 1 
a'a+ — 



(7) 



(8) 



One can then construct the Hilbert space as usual; the 
instantaneous ground state of the system at time t is 
defined by 



a\0 t >= 0, 



and an orthonormal basis of energy eigenvectors can be 
found as 



(at 1 )" 



(10) 



The operator equation (6) can be solved in terms of 
the constant operators at time to, i.e. a and aj, by 
introducing a Bogoligov transformation 



a = a ao + j3* aj , 
a) = (3 ag + a* aj. 



(11) 



The complex functions a and (3 obey the same differential 
equations as a and at given in (6), respectively. Eq. (7) 
requires 



M 2 -|/3| 2 = i, 



(12) 



and initially one should choose a(to) = 1, 0{to) = 0. 

The ground state of the system at time t , which is 
defined by 



a |0 >= 0, 



(13) 



will not remain to be the ground state at a later time. 
Rather, the expectation value of the number operator 
TV = at a in this state can be found as 



< N > = < 0o\N\0 o >= \(3\ 2 , 



(14) 



which shows that on the average this state contains |/3| 2 
quanta at time t. One can actually be more precise and 
determine |0q > exactly. Expanding it in the orthonor- 
mal basis vectors and using (11) and (13) it can be shown 
that up to an irrelevant phase 



|0o >= ^=cxp(^at a t)|0 t > 



(15) 



From this expression one can read of the probability Pm 
of finding an even number of quanta in |0q > 



In 



Pin — 
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(16) 



where the numerical prefactor can be recognized as the 
n'th order Taylor expansion coefficient of l/\/I — x. It 
is easy to see that Pi n is a decreasing function of n and 
for n » 1 one has 



1 



2n 



|/?l 



In 



^T|a| 2 "+ 1 ' 



(17) 



Although on the average |0o > contains |/3| 2 quanta, the 
most probable outcome of a measurement is the ground 
state 1 0t > with no quanta. 

To characterize the probability distribution one can 
calculate the deviation A7V of the number of quanta 



(AN) =< N > — < N > . 



(9) Using 



< N 2 >=< o |iV 2 |0o >= 2\a\ 2 \[3\ 2 



one finds 



AN r- a r- 

Ar - V2 |— I > V2. 



(18) 



(19) 



(20) 



For large average production, i.e. \/3\ » 1, the relative 
deviation is equal to y/2. In the opposite limit of small 
quanta creation, i.e. \(3\ <C 1, it is given by y/2/\/3\, 
which is much larger than unity. In any case, the main 
conclusion here is that the number of created quanta has 
large fluctuations about the mean value given by (14). 
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III. QUANTUM PARTICLE PRODUCTION: A 
WAVE PACKET ANALYSIS 

Consider a real scalar field X propagating in a cosmo- 
logical Robertson- Walker background 

ds 2 = -dt 2 + a 2 (dx 2 + dy 2 + dz 2 ), (21) 

which has the action 



s = -\Jv—g[(Vx) 2 + M 2 x 2 ]. 



(22) 



We assume that in addition to the scale factor a, the 
mass parameter M may also depend on time t. Defining 
a new field by 



X = a^ 2 X , 
the action up to surface terms becomes 



(23) 



X 2 



_ (M 2 - V - Ih)X 2 



where dot denotes time derivative and H = a/a is the 
Hubble parameter. One can expand the field in Fourier 
modes and introduce ladder operators as in the harmonic 
oscillator system discussed in the previous section 



X = 



1 



(2tt) 3 / 2 



d A k 



/2oj k 



ik.-x.-i j t iAj k dt _|_ ^ c 



where h.c. denotes hermitian conjugate, k is the comov- 
ing wave vector, k 2 = fc'fc- 7 ^ = k.k and 



,2 _ , ,2 , k2 9 rj2 3 ••• 



(24) 



uj 2 = M z + ^--H 2 - -H. 
k a 2 4 2 

The time dependence of the operators ak and at are im- 



posed to get the conjugate momentum II = X as 



it^fcOk -ik.x-i f* u k dt , , 

e J + h.c 

/2uj k 



which implies 



dt k = ^ e - 2i f"* dt a k . 
2uj k 



(25) 



In this case the equal time canonical commutation rela- 
tion [X, IT] = iS(x — x') is equivalent to 



[ok,4'l = S(k-k') 
and the Hamiltonian can be expressed as 



H 



I 



d s k 



a k a k + 2 



(26) 



(27) 



The instantaneous ground state is defined by 

a k \0 t >= 0, (28) 

and the Hilbcrt space can be build by acting with the 
creation operators on |0 t >. 

One thus rediscovers the well-known fact that the free 
field theory of x particles is nothing but an infinite col- 
lection of harmonic oscillators each of which is labeled 
by a comoving wave vector k and by the time depen- 
dent frequency (24). The only (minor) complication is 
that due to conservation of momentum ladder operators 
having wave vectors -l-k and — k are coupled by (25). 

To solve (25) one can introduce a Bogoligov transfor- 
mation 

flk = a k a k (t ) + f3la f _ k (to), 
al k = P k a k {to) + a* k al k (t ),. (29) 

where 

Wfe 2i f ' ui k dt n 

a k = ^ — e J Pk, 
2uj k 

Pk = ~ — e J a k . (30) 

ZUJ k 

Note that a k and j3 k depend only on the magnitude k 
and not on the direction. 

As in the harmonic oscillator problem, the ground state 
of the system at time to, 



flk(M|Oo >= 0, 



(31) 



becomes a multiparticle state at a later time. Similar to 
(15), it can be expressed in terms of the ground state 
|0 t > at time t as 



ioo >=n 



i 



exp(^aL k 4) |0 t >, (32) 
ak\ 2a k 



i.e. for each momentum mode labeled by k the probabil- 
ity distribution for the number of such particles contained 
in |Oo > is given by (16), where half of the particles have 
momentum k and the other half have — k. Note however 
that as it stands the above formula does not make sense 
and needs regularization. 

The mean number of modes contained in |0o > can be 
found by calculating the expectation value of the number 
operator 



^Vk = a k fl k> 



(33) 



which by (26) reads 

< A^ k >=< 0o|^Vk|0 o |/5 fc | 2 <5(0). (34) 

The infinity contained in the delta function can be inter- 
preted as 



*(0) = 



V 



(2tt) 3 



(35) 
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where V is the total comoving volume of the space. Defin- 
ing the comoving number density 



«k = 



V 



one obtains the familiar expression 
< n k >=< |n k |0 >= 



W 
(2tt) 3 



(36) 



(37) 



which is well defined. 

Although it looks like one manages to make sense of 
the infinity in (34), there is still an interpretation prob- 
lem. The momentum modes we are dealing with are com- 
pletely dislocalized and it is not meaningful to talk about 
the density of such excitations. Moreover, the calculated 
mean values should be explained in an ensemble picture 
but in our case there is only one realization of the sys- 
tem, i.e. in a given spacetime one can only make a single 
measurement on the mode labeled by k. As pointed out 
in the previous section, the result of a measurement may 
not be close to the average value since the deviation from 
the mean is large. Indeed the most probable outcome of 
a single measurement is the state with no particle. 

There is an alternative and physically more viable way 
of dealing with the above problem. The infinity appeared 
in (34) can be thought to arise due to the particle creation 
in an infinitely large space. A similar divergence also 
appears in the particle creation by black holes which was 
interpreted as the finite steady rate of emission for an 
infinite time as shown in [12] by analyzing the process 
using localized wave packets. Following [12] we introduce 
a new set of ladder operators as follows. Let e > be an 
arbitrary comoving momentum scale and introduce two 
vectors J = (ji, .72, jz) and n = (ni,n 2 ,n 3 ) with integer 
entries. Let 



ajn = j exp(-27rin.k/e) a k , 



(38) 



where we introduce a shorthand notation for a three di- 
mensional momentum integral 

f-Ol + l)e /-C72+1)* r(j 3 + l)e 

j j dkxdk 2 dk 3 . (39) 

/J Jj ie Jj 2 e Jj 3 e 

As we show below, oj n is the annihilation operator for a 
mode peaked around the comoving position 27rn/e with 
a spread 1/e in each direction. For small e, the mode can 
be thought to have a comoving momentum Je. The new 
operators obey 

[ajn, flK m ] = S JK S nm , (40) 

which justifies the identification of aj n and a\ n as ladder 
operators 

It is possible to invert (38) and express a k in terms of 
aj n as 



a k 



-Jj^ ^ exp(27un.k/e) a Jn ; jit < h < (j; + l)e. 

(41) 



Using this expression in the Fourier expansion of X one 
gets 



* = E 



J,n 



/jn(x) 



i J uj k dt 



+ h.c. 



(42) 



where the mode functions are given by 

jjn(x) = 1 J exp(-ik.(x - 27rn/e)). (43) 

Performing the integral it is easy to see that /j n is local- 
ized around x = 2nn/e with width 1/e. Moreover 



/ 



d 3 X / Jn / Km = Sj K 6„ 



(44) 



and 



E^( x )/jn(x')=<5(x-x'), (45) 

J,n 

which imply that these functions form a complete or- 
thonormal family. From the expansion (42) we see that 
the operator a\ n acting on the ground state produces a 
localized field quanta with the wavefunction /j n . 

Ignoring the zero-point energy, the Hamiltonian (27) 
can also be expressed in terms of the new ladder operators 



as 



#=-^E f Ld k exp(2Trik.(n- m)/e)a\ n aj m . (46) 

6 J,n,m"' J 



Therefore the Fock space states 



10* > 



(47) 



are not the eigenvectors of the Hamiltonian. Note that 
the definition of the ground state in (28) is equivalent to 
ajn|0 t >= 0. 

If u)k does not change significantly in the interval 
jit < ki < (ji + l)t for a given J, then it may be taken 
out of the integral in (46). In this case, the momentum 
integral yields e 3 i5 nm and the Hamiltonian becomes di- 
agonalized in this particular subspace of basis vectors. 
If this condition is satisfied for all J then all the Fock 
space states become (approximate) eigenvectors of the 
Hamiltonian. 

Mathematically speaking, Uk can be treated as a con- 
stant in (46) if t is small enough such that 



(48) 



For a massive field this condition is satisfied for all mo- 
menta if t <C m and for a massless field it is only obeyed 
in the momentum range k 3> t. Assuming (48) holds, 
the Hamiltonian turns into 



. J 



(49) 



■5 



where j 2 = J.J. Recalling that the sum over n can be 
viewed as a sum over cubic regions of side length 2ir/e, 
one can define the total number and the energy density 
operators as 



E W * °in fl Jr 



(27T) 



(50) 
(51) 



In analogy with (33), one can also introduce the number 
operator for a specific mode 



Aj n = aj n aj 



(52) 



Note that by (40) ./Vj n is an honest counting operator. 

Let us now analyze particle creation effects using these 
localized modes. The expectation value of JVj n in the 
ground state (31) can be calculated as 

<iV Jn >=<0 o |iV Jn |0o>=^|/3 fc | 2 . (53) 

Dividing with the volume of the region one gets the num- 
ber density as 



< n Jn >= 



1 



(2tt) 3 



(54) 



Although (53) and (54) look like (34) and (37), respec- 
tively, it is a lot easier to interpret them physically. First 
of all, since we are dealing with localized excitations it is 
meaningful to talk about the number density of modes. 
Secondly, unlike the previous expression (34), (53) is free 
of infinities. Finally, even in a single realization, i.e. in a 
given spacetime, it is possible to make sense of the mean 
values since one can make independent measurements in 
different spatial regions labeled by n and calculate the 
average of the outcomes. Note that due to translational 
invariance the mean values (53) and (54) do not depend 
on n. 

As discussed in the previous section, an important 
property of the particle creation effects in this setup is 
the existence of large deviations about the average val- 
ues. This clearly indicates a breakdown in the homogene- 
ity since the number of particles with a fixed quantum 
number J produced in different regions labeled by n will 
substantially be different. 

There are, however, two points one should be care- 
ful about before concluding that the homogeneity is de- 
stroyed. Firstly, the interactions between particles work 
for homogenization. Since the Hubble distance sets a 
limit for the range of local interactions in an expanding 
universe, these can be neglected if e is chosen to have this 
scale. Secondly, for small e the spacing between neigh- 
boring momentum levels is narrow and statistics might 
work to reduce the deviation by averaging over nearby 
levels. To take into account this point, one should calcu- 
late the deviations of the total number of modes n t or of 
the energy density p e , defined in (50) and (51). 



To determine those deviations we first note that 
straightforward manipulations yield 



< n £ >=< 0o|n e |0 >= ^3 J d 3 k\(3 k 



(55) 



and 



(An £ ) 2 =< |n 2 |0 > -(< |n e |0 >) 2 



(27 



For small e obeying \(3 k \ 3> ed\(3 k \/dk, the integrals can 
be approximated, e.g. by 



\Pk\ 



(56) 



Collecting terms and converting again sums into integrals 
one obtains the final result for the relative deviation 



An e 

<n £ > 



V2e 3 / 2 



Cfrf 3 fc|/3 fc | 2 KI 2 

Jd 3 k\f3 k \ 2 



1/2 



(57) 



Similarly the deviation in the energy density defined in 
(51) can be found as 



<p t > 6 f d 3 k\(3 k \ 2 oj k 



1/2 



(58) 



Let us recall that normalization requires |afc| 2 — \(3 k \ 2 = 1. 

It is possible to understand the origin of e dependence 
in (57) and (58) as follows. In a homogeneous background 
the number of modes in a volume V is proportional to 
V and statistically the relative deviation of a total quan- 
tity is expected to decrease with the square root of the 
number of modes, i.e. like 1/y/V. In our case V ~ 1/e 3 
which explains e dependence in the above formulas. The 
deviation in each excitation labeled by J can be much 
larger by (20), and thus (57) and (58) set a lower bound 
in the degree of inhomogcneities on the comoving scale e 
which can only be saturated if local interactions are very 
efficient on that scale. 

To determine the power spectrum P(k) corresponding 
to (58), a Gaussian or top-hat window function can be 
introduced to probe the scale k ~ e (see, e.g., [13]). A 
simple calculation shows that (58) has a white noise spec- 
trum similar to sub-horizon thermal fluctuations studied 
in [14]. 



IV. APPLICATION TO REHEATING 

In this section, we consider the particle creation ef- 
fects during reheating in a scalar field driven inflationary 
model. The evolution of the background fields, i.e. the 
metric and the inflaton, is governed by the Einstein and 
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the scalar field equations 
8tt 



H 2 



3M| 



U 2 + V(<f>) 



■ dV 

<t> + 3H(j) + — = 0, 



(59) 



where V{<f>) is the scalar potential. In some classes of 
inflationary models, e.g. in chaotic inflation, the scalar 
oscillates about the minimum of the potential during re- 
heating, so one can take 



V 



1 



-rn 



2 ±2 



(60) 



where m is the inflaton mass. In a chaotic inflationary 
scenario V can be assumed to have exactly this form. 
Due to the expansion of the universe the amplitude of 
the oscillations gradually decreases in time so one can 
write the solution for the inflaton as 



(f) = <&(t) sin(mi). 



(61) 



Assuming $ <§; m<I>, the Einstein and the scalar equa- 
tions become 



which can be solved as 



a = ao 



to 



2/3 



3irmt 



(62) 



(63) 



As it is well known, the stress-energy-momentum ten- 
sor corresponding to the coherent inflaton oscillations is 
equivalent to the one for the pressureless dust. In this 
case, the combination 9H 2 /4 + 3H/2 which appears in 
(24) is equal to zero. 

We first consider particle creation effects during pre- 
heating in a chaotic inflationary scenario with a quadratic 
potential (60). This has been studied in detail both an- 
alytically and numerically in [10] and we mainly use the 
findings of that paper below. In such a model inflation 
occurs when <f> > M p and it ends when the field decreases 
below 4> ~ Mp/2. It turns out that even after a single 
oscillation the amplitude drops enormously and the solu- 
tion (63) becomes a very good approximation. Following 
[10], we take to = ir/2m as the time for the beginning of 
preheating and thus initially 



$ 



2M 



v 



7rv37r 



— ' H ° ~ 3^ (64) 



A realistic value for the inflaton mass is m = 10 6 M, 



v 



We assume that inflaton is coupled to a light boson \ 
with an interaction term 



1 



c mt = --g 4> x , 



(65) 



where g is a dimensionless coupling constant usually as- 
sumed to be small. In this case, the results of the previous 
section can be applied by fixing the time varying mass 
in (24) as M 2 = g 2 <& 2 sin 2 (mt). As shown in [10], for 
9 3> 10~ 6 one finds a decay due to broad parametric res- 
onance and the occupation numbers of the modes which 
have comoving momenta less than fc* are exponentially 
growing, where 



— = \/gm$ Q . 
a 



(66) 



The particle creation effects in this model are charac- 
terized by an effective index \i k such that 



k = e^ mt . 



The index can be estimated as 



1 



P-k — M — oMfe {km) {k ~ fcm) 



(67) 



(68) 



where the second derivative of fik can be approximated by 
Mfe(fc m ) = 2fi/Ak. Here k m and Ak, which are the max- 
imum and the width of the first and the most important 
resonance band, are nearly equal to fc*/2. Depending on 
g the average index \x varies between 0.1 and 0.2 and for 
numerical estimations one can take /i ~ 0.13 [10]. 

Assume that the resonance ends at time t\ after N 
oscillations, where mt\ — 2irN. From (62) and (63) the 
Hubble constant at the end of resonance can be found as 
Hi = m/{3nN). In calculating the relative deviations, 
we choose e to be the comoving Hubble scale at the end 
of preheating, i.e. 



aiHi. 



(69) 



It is easy to see that in the broad resonance regime e <C fc* 
and the condition (48) is satisfied. 

Using (68), the integrals in (57) and (58) can be eval- 
uated by the steepest decent method [10]. Nothing that 
in this regime a k — Pk, one can straightforwardly find 



_Ap ± _ 

< Pe> 



An e 

<n f > 



2 9/4 
337^2 



/_^\V 4 fg^o 



-3/4 



For small g, the resonance may end before back reaction 
and rescattering effects become important. The number 
of oscillations for the first stage, where these effects can 
be ignored, are approximately given by N ~ gM p /(6imi) 
[10]. Usually, the time for the second stage is much more 
shorter so one can take N to have this value in estimating 
the deviation (70). Setting 
finally obtains 



.13 and m = 10~ 6 M p one 



Ap e 



An e 



2 x 10~ 7 



< p e > < n e > 



(71) 



For g <C 3 x 10~ 4 resonance ends in the first stage and 
the relative deviation is much larger than 10~ 3 . For big- 
ger g, one should consider back reaction and rescattering 
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effects and the deviation (71) gives the degree of inho- 
mogeneity in the beginning of this stage. Note that even 
for g = 10 -2 the relative deviation has the same order 
of magnitude with the fluctuations in the CMB temper- 
ature. 

As pointed out above, (70) actually gives a lower bound 
since the deviations in the numbers of individual excita- 
tions can be much larger. By choosing e to be the Hubble 
scale at the end of preheating we make sure that local in- 
teractions cannot restore homogeneity. 

Depending on the details of the theory, the reheating 
process may continue after preheating but the decay of 
the inflaton ceases to exist in the parametric resonance 
channel. In this period, the well known perturbation 
theory of reheating is applicable to study the particle 
creation effects. We now continue with the determination 
of relative deviations in this regime in a different model. 

In the following we assume a trilinear coupling, which 
may arise after spontaneous symmetry breaking, 



£ir, 



1 



(72) 



where a is a constant with mass dimension. As shown 
in [15], the preheating picture completely changes when 
both interections (65) and (72) present in the Lagrangian. 
Therefore, the perturbative decay due to (72) should be 
considered on its own, i.e. it is not to be preceded by 
the preheating considered above. Our aim here is to de- 
termine fluctuations in a well-known perturbative decay 
scenario based on (72). 

The particle creation process via trilinear interaction 
was outlined in [9]. The perturbation theory is applica- 
ble when aQ/m 2 <C 1. Note that the evolution of the 
background fields is still given by (63), however the am- 
plitude <E> is now small compared to its magnitude in the 
parametric resonance regime. The frequency (24) corre- 
sponding to the interaction (72) is given by 



^1 



k 2 

+ <7$ sin(mi). 



(73) 



Due to the expansion equivalent to matter domination, 
H dependent terms in (24) cancel each other. 

For |/3fe| <C 1, an iterative solution to (30) is given by 



Pk 



dt — cxp -2i / u>k(t')dt' , (74) 



where to and t\ denote the beginning and ending of the 
decay, respectively. This integral can be evaluated using 
the stationary phase method [16]. It is easy to see that 
uJk/uk term in (74) consists of a non-oscillatory ignorable 
piece together with a term proportional to cos(mi) . Since 
the phase in (74) is negative definite, only e +rmt part 
contributes to the integral. As a result one finds 



Pk 



m at, 



(75) 



where the phase function is 



f(t) = J [-2uj k (t')+m]dt'. 



(76) 



For a given k, the main contribution to the integral comes 
from an interval near t* fixed by u>k{t*) = rn/2, which, 
in the perturbative regime cr$/m 2 <C 1, implies 



k m 
a* 2 ' 



(77) 



This corresponds to the decay of the inflaton at time t* 
to two x particles with comoving momentum k. 

To apply the stationary phase approximation we first 
note that /(i*) = mH{t*) — 2<r$(f*) cos(mt*). From 
(62) and (63), the oscillating factor can be ignored for 
m 2 /M p > a which yields 



Pk^ 



2m 



I exp[imH(t*)(t-t Sf ) 2 /2] dt. 

J t 



Since m 3> H, the stationary phase approximation is 
valid and the limits of the integral can be extended to 
infinity. This up to an irrelevant phase gives 



a Tr<I>(t*)M p 



Pk — 

m z V 2 

To determine $(£*), we note that (63) implies 

3/2 -- -m\3/2 



2k) 



(78) 



(79) 



where in the last line we use (77) to determine a* in terms 
of k. This gives 



|2 7R7 2 M P $ fa \ 3 / 2 

\pk\ - 



2m 5 / 2 



\2k) ' 



(80) 



where the comoving momentum k should be in the decay 
range 



aom , aim 
< k < —— 



(81) 



to 



with a\ being the scale factor at time t\. In the perturba- 
tive regime cr$/m 2 <C 1 (and for m 2 /M p > a), \Pk\ <^ 1 
which justifies both the solution (74) and the station- 
ary phase approximation used in the calculation of the 
integral. 

The whole process ends when the energy density p x of 
X particles cathes up the energy density of the inflaton, 
i.e. 

where $i is the amplitude of the oscillations at the end 
of the decay. Using (80) and recalling that the integral 
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in the momentum space is defined in the spherical region 
(81), we find 



Oi 



5/2' 



(83) 



For ai > a , (82) gives $i - a 2 M p /m 2 . 

It is interesting to compare (83) with the equation 
(8.32) of Kolb and Turner [13], which gives the evolu- 
tion of the radiation energy density extracted from the 
inflaton by a constant decay rate T. It is easy to see 
that both equations agree for T <~ <r 2 /m, which up to a 
numerical factor is the total decay rate corresponding to 
the interaction (72). 

Using the value of <&i determined above, the Hubble 
constant Hi at the end of the decay can be found from 
(62) as Hi ~ a 2 jm. Since our aim is to calculate the 
deviation in the number density of particles produced in 
a Hubble volume at the end of reheating we choose 



aiHi ~ ai — . 

m 



(84) 



For m 3> <J and k in the decay range (81), e <C k and 
thus (57) and (58) is valid. Using (80) and noting that 
at the end of the decay uJk — k/ai the deviations can be 
estimated as 



Ap e 

< Pe > 



An e 

<n e > 



M„ 



(85) 



Assuming m = 10~ 6 M P and a = 10~ 7 m (to satisfy 
the requirement m 2 /M p > a imposed by the stationary 
phase approximation), one finds that the relative devia- 
tion is of the order of 10~ 13 , which is very small compared 
to the deviation found after preheating. 

The reheating temperature in this model can be esti- 
mated from the energy density at the end of the decay 
(see e.g. [1, 9]) 



m 2 $ 2 



(J Mr, 



m 



(86) 



which gives a very low reheating temperature Tr ~ 10 9 
GeV. Therefore the decay actually happens rather slowly 
giving a large spectrum of created particles to be av- 
eraged over and a very small e. This explains both 
the smallness of the deviation and the contrast with 
the quantum mechanical result (20), which gives an ex- 
tremely large deviation in perturbation theory. Note that 
the reheating temperature is independent of the initial 
value of the amplitude <&o, as it should be. 

To determine the deviation for larger reheating tem- 
peratures, we first note that in the perturbative regime 
|cKfe| ~ 1 and (57) becomes 



< n e > 



V2e 3 / 2 



(Jd^k\p k \ 2 y/ 2 - 



(87) 




The integral in (87) can approximately be estimated in 
terms of the reheating temperature as T^/m, where we 
set oi = l. Similarly the Hubble scale at the end of 
reheating is H ~ T^/M p = e. Therefore, the relative 
deviation in the number density of particles produced in 
a Hubble volume is 



(88) 



For T R = 10 15 GeV and m = 10~ 6 M p the relative de- 
viation becomes 10~ 7 . Thus the deviation appeared in 
perturbation theory is gcncrically smaller than the one 
encountered in the parametric resonance regime. This 
result can be explained by the fact that preheating oc- 
curs suddenly in a very short time (i.e. explosively) and 
in a comparatively narrow band as measured by the Hub- 
ble scale. As a result the number of the modes to be 
averaged over is small as compared to the perturbative 
regime, which lessens the statistics and gives a larger to- 
tal deviation. 



V. CONCLUSIONS 

In this paper, we study the particle creation effects in 
the context of quantum fields in time dependent external 
backgrounds. First considering the quantum mechani- 
cal harmonic oscillator with a time dependent frequency, 
which is the prototype of the field theory problem, we 
observe that the number of created quanta largely devi- 
ates about its mean value. Armed with this observation, 
we use a complete orthonormal family of localized wave 
packets to calculate the deviations in the field theory side. 

The wave packets that we introduce following [12] are 
defined with a fixed momentum scale e. Roughly speak- 
ing, they carry discrete momentum values which are split 
up by e and are localized in a region of size 1 /e. For small 
enough e, they also become eigenvectors of the Hamilto- 
nian. Therefore, these wave packets form a legitimate 
basis in the Hilbert space as good as the Fourier modes. 
Indeed, as discussed in section III, it is a lot easier to 
interpret them physically as long as the measurements 
are concerned. 

To understand e dependence of the relative deviations 
assume that particles are produced around momentum k 
in a band having a width Ak with equal average produc- 
tion. For each wave packet in the band, the deviation 
in the number of created quanta is given by the quan- 
tum mechanical result (20). There are however nearly 
k 2 Ak/e 3 number of different wave packets in this band. 
Therefore, the relative deviation of the total number 
of created particles should decrease like e 3 / 2 / '(k 2 Ak) 1 / 2 . 
The deviation found in preheating (70) can easily be re- 
derived by noting that Afc ~ k <~ y/gm&o and e <~ to. 
Since the decay is very quick the dependence on the ex- 
pansion of the universe is very weak in this case. How- 
ever, in perturbation theory the decay is very slow and 
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this naive estimate should be refined. Note also that de- 
creasing e is equivalent to observing a larger space and 
thus it is natural to get a small deviation. In an experi- 
mental setup, the e parameter should be set by details of 
the measurement. On the other hand, in an expanding 
universe the Hubble distance gives a natural scale. 

Applying this general framework to inflation we find 
that there emerges density inhomogeneities on Hubble 
length scales at the end of reheating. Depending on the 
details of the decay, the relative order of these inhomo- 
geneities can be as large as 10~ 3 in preheating and 1CP 7 
in perturbation theory. Note that local interactions are 
incapable of restoring homogeneity since they cannot op- 
erate on distances larger than the Hubble radius. More- 
over, these estimates actually give a lower bound since 
the deviation in the number of each particular mode can 
be larger. On the other hand, it is much safer to esti- 
mate the deviation in the total number of created parti- 
cles since for small e the nearby momentum levels become 
physically indistinguishable. 

It would be interesting to sec whether these density 



inhomogeneities affect cosmological events following re- 
heating like baryogencsis. It is also important to deter- 
mine if they switch to the non-linear regime and grow in 
time due to gravitational collapse (note that since they 
are different than the usual Fourier modes, the standart 
arguments cannot be applied here) . If this does not hap- 
pen, interactions can gradually restore homogeneity as 
the Hubble radius grows, but the time required for full 
equilibrium might considerably be large which would re- 
duce the reheating temperature. In any case, due to the 
quantum mechanical nature of the particle creation pro- 
cess in reheating, we find that the hot stage does not start 
with the perfect homogeneity attained after inflation. 
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